arXiv: 1501.04512v2 [math-ph] 13 Feb 2016 


From continuum mechanics to SPH particle systems and back: 
Systematic derivation and convergence 

Joep H.M. Evers* Iason A. ZisisT Bas J. van der Linden! Manh Hong Duong§ 


Abstract 

In this paper, we derive from the principle of least action the equation of motion for a 
continuous medium with regularized density field in the context of measures. The eventual 
equation of motion depends on the order in which regularization and the principle of least 
action are applied. We obtain two different equations, whose discrete counterparts coin¬ 
cide with the scheme used traditionally in the Smoothed Particle Hydrodynamics (SPH) 
numerical method (e.g. j22|h and with the equation treated by Di Lisio et al. in [J], 
respectively. Additionally, we prove the convergence in the Wasserstein distance of the 
corresponding measure-valued evolutions, moreover providing the order of convergence 
of the SPH method. The convergence holds for a general class of force fields, including 
external and internal conservative forces, friction and non-local interactions. The proof 
of convergence is illustrated numerically by means of one and two-dimensional examples. 

Keywords: Smoothed Particle Hydrodynamics, principle of least action, Wasserstein 
distance, measure-valued equations, convergence rate 
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1 Introduction 

The Smoothed Particle Hydrodynamics (SPH) numerical method was initially introduced to 
solve the equations of astrophysical flows. In the course of time it found application to equa¬ 
tions describing a plethora of physical processes (for its diverse applications, see [23]). These 
processes predominantly involve continua and the equations refer to systems with infinite 
degrees of freedom. The central idea of the SPH method is to set up a relation between 
the continuum and a particle system, in which the continuum is loosely considered to be the 
limit case in which the number of particles tends to infinity. Here, a ‘particle’ should not be 
interpreted as a physical object of any scale (like an atom, molecule or grain) but rather as a 
numerical entity attributed with mass, position, velocity and other properties of the medium 
it represents. 
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It is well-established that the classical SPH scheme can be derived formally by applying 
the principle of least action to the particle system, where the SPH density approximation 
acts as a constraint; see e.g. mmm ■ The importance of the particle system’s Lagrangian 
function was already recognized in the first articles describing SPH; cf. [16J. A subtlety lies 
in the fact that in the derivation of the SPH equations, the action of the particle system is 
minimized rather than the action of the continuum. The minimization of the action at the 
continuum level and the subsequent discretization of the motion equation in terms of particles 
do not necessarily yield the same equation (at the discrete level). 


The main achievement of this paper is twofold: 

• We introduce a systematic procedure for deriving measure-valued and particle formu¬ 
lations of continuum mechanics equations. We obtain two different schemes depending 
on the stage at which a regularization of the density is introduced. See Section [2] 

• We prove the convergence of both schemes using the Wasserstein distance on the space 
of probability measures; cf. Section [3j 

We now describe the two parts of our paper in more detail. 


In the first part (Section [ 2 ]), we aim at clarifying the exact difference in outcome between 
minimizing the action of the particle system and minimizing the action at the continuum 
level. To achieve this, we introduce a systematic procedure consisting of the following three 
steps: 

A formulation in terms of measures and, simultaneously, the regularizatioi^]of the density; 

B introduction of a particle formulation; 


C application of the principle of least action. 


These three steps are introduced in more detail in Section 2.2 It turns out that the order in 
which these steps are executed determines what the resulting equation is. To be more precise, 
the classical SPH scheme (as described e.g. in [223 ) is obtained, whenever the regularization 
of the density takes place before applying the principle of least action. That is, whenever the 

PTP1 


steps are executed in the order ABC 


or 


Both procedures are presented here; see 
Sections |2.3| and 2.4 If we apply the principle of least action (to the action at the continuum 
mechanics level) before regularizing the density then we obtain a scheme that appears in Di 
Lisio et al. [3 and in the recent paper [6]. However, this variant of the scheme is studied far 
less in literature. The procedure to obtain this scheme follows the order [C]|Ap 3| Its distinct 
characteristic is that it requires the gradient of the pressure field to be expressed analyti¬ 
cally, while the pressure itself does not appear in the numerical scheme, in contrast to the 
commonly used SPH schemes. We emphasize that although both schemes arrive from the 
principle of least action, the latter can also be derived directly from Newtonian mechanics 
and introduction of the density regularization. The details of our rational derivation and 


lr rhe regularization of equations is an old concept, introduced by Friedrichs in 1944 m- Additionally, 
notice that the regularization kernels used in SPH are a special subclass of the modifier functions used by 
Friedrichs the SPH kernels are symmetric positive mollihers. 
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the mutual relation between the two schemes have, to our knowledge, not been described in 
literature before. 


In any case, regularization of the density practically means that the original problem is de¬ 
liberately turned into a regularized one, which is afterwards solved by means of some variant 
of an SPH scheme. Hence, by choosing SPH as the solution method one is automatically 
bound to studying a different problem than the original one at the continuum level. Thus, 
two questions naturally arise: 

• Does the solution of the regularized problem converge to the solution of the original 
problem? 

• Does the particle solution of the regularized problem converge (in a certain sense) to 
the solution of the regularized continuum problem? 

The former is out of the scope of the present study (|S] has dealt with it), while the latter is 
the topic of the current work’s second part. 


Measure theory provides a framework to study the limiting behaviour as the number of parti¬ 
cles goes to infinity (cf. also e.g. [21]). Both the particle system, and the limiting continuum 
setting can be formulated in terms of measures. Hence, a distance between measures is a 
natural tool to characterize convergence; in this work we take the Wasserstein distance on 
the space of probability measures. This particular distance has the advantage that it can be 
formulated as the infimum over a set of joint representations (more details follow in Defini¬ 
tions 3.2 and 3.3). This is convenient, since one can thus obtain an upper bound (needed to 
prove convergence!) by choosing any admissible joint representation. See also [30], Chapter 
6, for more discussion. 


We prove the convergence of measure-valued solutions, as the initial measure is approxi¬ 
mated; cf. Section [3] The line of arguments is similar to the one followed by Di Lisio et al. in 
who first employed measures in combination with the Wasserstein distance to prove the 
convergence of the SPH method, but the result obtained in the present work is more general. 
It should be mentioned that in the earlier work |24| convergence of the empirical measure 
representing the particle system was proven, but using a different technique. Moreover, the 
only forces considered were mutual interactions between particles. Other approaches to ob¬ 
tain convergence are given e.g. by [2] using maximum local entropy estimates, )26| employing 
estimates for the truncation error, and EH H8]. 

Nonetheless, the scheme treated in |7] is not the aforementioned traditional scheme. Our 
proof applies both to traditional SPH and to the scheme covered by [7]. Moreover we allow 
for a much more general class of force fields, including external and internal conservative 
forces, as well as friction and non-local interactions. 

Previous work in the framework of measures by the authors of the current paper can be found 
in eh. where apart from the aforementioned force terms also random noise is treated. In 
E3H3 measure-valued evolutions are treated in the scope of equations of motion that are 
first-order in time. The link between first-order and second-order models is discussed in [12]. 


The theoretical result of this paper regarding the order of convergence is supported numeri¬ 
cally in Section 3.6 for one and two-dimensional illustrative examples, which involve different 
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force fields. 


In Section [4] concluding remarks are given about Sections [2] and [3j Also, some attention 
is given to possible future research directions. 


2 Systematic derivation of the equations of motion 


In this section we derive equations of motion from Hamilton’s principle of least action, which 
involves the Lagrangian function posed in a continuum mechanics setting. We describe an 
explicit ‘recipe’, hence avoiding the need to introduce approximations in an ad hoc manner. 
This ‘recipe’ consists of three building blocks (coined [A| [B| and [Cj see Section 2.2). The order 
in which these blocks are executed, influences the final outcome. As such, the systematic 
procedure we describe here also shows exactly how different formulations/schemes arise from 
the same basic principles. 


2.1 Derivation of the action in a continuous setting 

Assume that for fixed time t a mass density pt on a spatial domain Lit is given. We define 
the Lagrangian density of our system as 

£{p-,V,u) '■= Qm 2 -e(p{y),y)^j p(y), (2.1) 

where y and u are independent Eulerian coordinates, and e denotes the internal energy density. 
To obtain the Lagrangian L, we integrate C over the spatial domain Lit.: 


L(t) := f C(p t ;y,u)dy. (2.2) 

fit 

For this integration to make sense, we assume now that u is actually a velocity field, defined 
as a function of t and y: u := u(t,y). Let there be a coordinate transform such that 
Lit = < ht(Ho) for some initial domain Do- We call the family of transformations ( < hf)t^o a 
motion mapping and transform the integral above according to y = <hj(x) with x G LIq: 


m(t) = 


fio 


-\u(t, $t(x’))| 2 - e(pt($t(x)), $t(x)) ) Pt($t{x)) \J$ t (x)\ dx. (2.3) 


The functional dependence of L on the motion mapping is indicated by explicitly including 
in square brackets. The expression | Jd*t| denotes the determinant of the Jacobian matrix 
of the transformation, consisting of the derivatives of the components of ( 3L with respect to 
the components of x. Now we assume that the density pt relates to the density po defined 
on the original domain LIq by the same transformation which is mathematically described 


by a push-forward , pt = &tffPo (cf. Definition 3.1). In particular, the densities relate in the 
following way (see e.g. 0, P- 90): 


po(x) = pt($t(x)) |J$t(s)| 


(2.4) 
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Combined, (2.3) and (2.4) yield 


m(t) = 


f!o 


| u(t, 4h(x))| 2 - e(pt($t(x)), 4>i(x)) p 0 (x) dx. 


(2.5) 


In the above we fixed t, but obviously all arguments can be repeated for every t in some 
interval [0, T]. In particular, we are interested in those motion mappings that are continuous 
and differentiable in time, and we wish to obtain their equation of motion. The introduction 


of the motion mapping (4>i) te [ 0 ,T] has taken us from pure Eulerian coordinates in (2.2) to¬ 
wards Lagrangian (material) coordinates in (2.3). The crucial and final step to complete this 


procedure is now to specify what the velocity field u is. In order to remain consistent with 
the motion mapping we introduced, we postulate the relation: 


u(t, $t(x)) = <&t(x) for all x G JIq- 


( 2 . 6 ) 


The velocity u(t, 4h(x)) is the velocity at time t of a material point that started in x at time 
0, and — in words - ( |2.6[ ) means that this velocity is equal to the time derivative at time t of 
the position 4h(x) of that particular material point. By connecting the Eulerian velocity u 
to the Lagrangian velocity 4>t, we obtain the Lagrangian functional 


L[4>](t) = J ^|4> t (.T))| 2 - e(p t ($ t (x)), 4h(x))^ po(x) dx. 

fio 

We define the action of our system by 



L[4>](f) dt. 


(2.7) 


( 2 . 8 ) 


2.2 Three procedures 


The aim of this part of our paper is to derive equations of motion from the action (2.8), by 
means of the Euler-Lagrange equations (we will see that these appear in different shapes). 
Moreover, we wish to derive these equations of motion for a particle system, which naturally 
induces a numerical scheme. A methodological way to go from the continuum (Section |2.1| ) 
to a particle system, is via a measure-valued formulation. Our motivation to do so is the 
fact that we need a framework that incorporates the ‘real physics’, i.c. the density pt, and 
an approximating particle system to establish the convergence of the particle scheme to the 
continuum. 


To get the transition from the continuous action (2.8) to equations of motion for the par¬ 
ticle positions, three steps are necessary: 


A introduction of measures: replace pt(x)dx by pt(dx) and, wherever necessary, approxi¬ 
mate pt by some pt that depends on 

B substituting for p t a discrete measure p,™ = 

C Derive the Euler-Lagrange equations (either classically or in variational sense). 
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The steps are here described in a somewhat simplistic and unprecise way; their true meaning 
will become clear in Sections 245, 2 A and |2.5| Step [A] takes us to a regularized version of the 


problem, which is a problem different from the original one. Step B cannot happen before 


A, but we have the freedom to choose the further ordering. This gives rise to three different 


derivations: 

this procedure discretizes the Lagrangian and derives the corresponding equations of 


motion afterwards; see Section 2.3 


Kirn this procedure derives the equations of motion from the measure-valued Lagrangian 


and discretizes these equations afterwards; see Section 2.4 


\nm this procedure derives the equations of motion from the continuum Lagrangian, writes 
them in measure-valued form and discretizes afterwards; see Section 12.51 


Procedures 000 and |A||C|p| eventually yield the same particle scheme. This is the scheme 
traditionally used in the SPH community (cf. |22j). Procedure C A^ is the one that yields 
the equations used in |7J and [6]. 
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2.3 Equations of motion via the route AjBflC 
2.3.1 Step [A] 


In Section |2.1| we introduced (for each t ) the density pt as the push-forward of the initial 
density po under the mapping <E. In this section we lift the evolution of pt to the space of 
(time-dependent) measures. Let po and pt be the measures associated to the densities po and 
p t . Hence, pt = &t#P o- hi (2.7)—(2.8) we can substitute po(x)dx by po(dx). Afterwards, 
there is one more aspect that we need to ‘repair’ before we are completely in a measure 
formulation. The internal energy density e depends on pt itself, via pointwise evaluation at 
<L t (x). An approximation of pt is needed to obtain a general expression that is even well- 
defined for measures that have no density (w.r.t. the Lebesgue measure). We propose to 
introduce a regularization via convolution 


Pt(0 ■■= (W h * Pt) (0 = fW h (t - y ) pt(dy), ( 2 . 9 ) 

n t 

for all £ £ M d . Here, the smoothing function Wy L is nonnegative and even (so that it obtains 
an odd gradient, an effect which is used later in the derivation of the equations), h is a small 
parameter, and Wy t —*■ So in the narrow topology as h —> 0 (i.e. tested against bounded contin¬ 
uous functions). A typical example is the Gaussian with zero mean and variance h 2 / 2 . If pt 
has a density pt then the convergence pt —> pt. holds in some sense and under certain mathe¬ 
matical conditions. E.g. if p t is continuous and bounded, then by definition of Wy —*• So, ptiO 
converges to pt(£) for all £. In any case, the convolution regularizes the solution, introducing 
an artificial ‘density’ pt, such that pointwise evaluation and the gradient are defined even 
when pt does not exist or is not differentiable. Note that, pt also depends on h, but in this 
work we do not consider the limit h —> 0, therefore for simplicity of notation, we leave out h 
in pt . However, we stick to the subscript h in Wh in agreement with the common notation in 
SPH literature. 


Note that pt can also be written as 

Pt(0 = [W h (d - $ t (z)) po(dx), (2.10) 

by definition of the push-forward. Hence, we should keep in mind that p has either a func¬ 
tional dependence on pt, or an extra dependence on $i(-) (depending on which formulation 
we choose), but we do not write this dependence explicitly. 

In e, we substitute pt for pt in the sequel and redefine the Lagrangian (in a measure- 
formulation) such that the action becomes 

T T 

Ql^tOOl 2 -e(/5t($i(z)),$t(z))^ po(dx)dt. (2.11) 

o o n 0 

The new, generalized formulation in terms of measures allows us to consider more types of 
solutions, simply by allowing for more general initial conditions. This is exactly what we 
exploit in the following step via a particle approximation. 


sm = Jm(t)dt=J J 
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2.3.2 Step [B] 

In this step, we substitute for /iq a discrete measure of the form rn id Xt 0 - Under 

push-forward, the measure remains a discrete measure with positions of the Diracs {xj(f)} 
evolving under the motion mapping: Xi(t) = o)- We emphasize that the equation for 4^ 

is yet unknown and is to be derived in the next step. 

The Lagrangian takes the form 


L [®]{t) =i>2 mi Ql^WI 2 > (2-!2) 

with 

n 

Pt(xi(t )) = Y,mjW h (xi(t) -Xj(t)). (2.13) 

3 =i 

In the literature of SPH, particles of the same mass are employed for the modeling of the flow of 
a single fluid. In that case, the term rrii corresponds to 1/n. On the other hand, multiphase 
media of piecewise continuous mass density can be modeled with the use of particles of 
different masses [221EU. For that reason, we adopt the general case of (in principle) unequal 
masses m*. 


2.3.3 Step [C] 

The equations of motion are obtained via the ‘classical’ Euler-Lagrange equations, see (1.57) 
in m, applied to the Lagrangian 


i —1 



W h {yi - yj ), yi 


D=i 


(2.14) 


cf. (2.12). In the presence of nonconservative forces (cf. p. 23 in |17| ) the corresponding 
Euler-Lagrange equations are 


dt 


V Uk L 






= m k q[p%]($ t (x k , o), $ t (x k>0 )), (2-15) 


for each k G {l,...,n}, where q is the force density (per unit mass) of nonconservative 
forces. The functional dependence in square brackets denotes that q incorporates a non¬ 
local interaction term. More details will follow later; cf. (2.34). The subscript “(yi,Ui) = 
(^(a^o), $ t (x ii0 ))” should be read as performing this substitution for all* G {1,..., n}. 
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After calculating the derivatives V Ufc and X7 yk in (2.15), we obtain 


m k — $t(xk, o) = 

de I n 1 n 

~ mk !f mj Wh ^^ Xk d) ~ ^( x i,o)), ®t{x k , o) ) ^2 m j VW h ($t(x k , o) - ®t{xj,o)) 
\j =i ) j =i 

n q I n \ 

+ y~! m J w h{$t(xi, o) - ^(x j>0 )),^Ko) 1 m fc Vlff l ($t(2:j i o) - 4> t (x fcj0 )) 

i=l V 7=1 / 


v-7. 

n 


-m k V y e Wft(^t(xfe, 0 ) - $t(®j,o)),$t(®fc,o) I + m k q[pt]($ t (x ki0 ), $t(xk,o))- 

(2.16) 

We denote by V y e the gradient of e only in the explicit spatial coordinate; that is, the second 
variable of e. We divide all terms by m k (which is nonzero without loss of generality). If in 
the second line we take del dp inside the sum and we use in the third line that Vlf/, is an 
odd function, then the corresponding terms in (2.16) can be combined, and we obtain 


&t(xk, o) = 

n 

- 22 mi ^ W h{$t(x k ,o) ~ $t(Xi, 0 )) 

2=1 

- V y e (pt($t(x kfi )), $t(x kt0 )) + q[p2]($ t (x k , o), $t(zfc,o)), (2.17) 
for each k E {1,..., n}. For brevity of notation, we use p again in the argument of e. 


de de 

— (pt(&t(x k , o)), ®t(x k , o)) + 7^ (pt(®t(Xi,o)),®t(.Xifl)) 


2.4 Equations of motion via the route |A|P[B 
2.4.1 Step [A] 


This step is exactly the same as in Section 2.3.1 


2.4.2 Step [C] 


We start from the action given in (2.8). Instead of using the classical Euler-Lagrange equa¬ 
tions, we employ here a generalized form of the principle of least action (see p. 127 and Section 
4.4 of [5]): 

S'[$](tt) = -<?[$](*), (2.18) 


for all test functions T E C£°((0, T); C'(?°(f2o;M d )). Here, S ,, [<h](4') denotes the variational 
derivative of S in the direction of T, and <5[4 > ]('k) is the work done along T. It is defined as 


T 

Q[$]W := J jq[pt]($t(x),$t(.x)) ■ ^t(x) p 0 {dx) dt, (2.19) 

0 Qo 
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where q is the force density as in (2.15). For S' we have: 




£=0 


Note that 


de 


e / W h (M x ) + £^ t (x) - $ t (y) ~ £^t{y)) Po(dy), M x ) + eM x ) 


£=0 


de 

dp 


W h (M x ) - My) Po(dy),M x ) / VWft($ t (i) - ^i(y)) • (^t(x) - %(y)) p 0 (dy) 


Uo 


+ v„e | y W fc ($t(x) - My) Po(dy), Mx) ■ Mx)- (2.20) 


To avoid lengthy notation, we denote here by e / [<h](^')(x) the expression in (2.20). The 
variational derivative of S can be expressed as: 


( 2 . 21 ) 


S'[$\ (tt) = / / (Mx) ■ Mx) - e / [^](^)(a;)J p 0 (dx) dt 
~Mx) ■ ’i’t(x) - e'[$](T)(.x)) p 0 (dx) dt, 


o n 0 

T 


0 fto 

where the last step follows from integration by parts with respect to the time variable. The 
boundary terms disappear because \h has compact support within (0, T). 


We rewrite the part involving 'Ff(y) in (2.20) as follows: 

T 

j j ~Wp (^(^i^))’^^)) J VWhiMx) - My)) ■ ^t(y) Po{dy) p 0 (dx) dt 

0 fig Qq 

T 

= IIJy p (pt(®t(y)),My)) v w h (M x ) - Mv))Mdy) ■ M x )po(dx)dt, (2.22) 

0 f2o ^0 

by subsequently interchanging the order of integration, using that the function VW is odd, 


and replacing x by y and vice versa. A combination of (2.18), (2.19), (2.20), (2.21) and (2.22) 
yields for 5 7 [d>](T) + Q[$](\li) an integral of the form 


...] • Ti(x) po(dx) dt, 


(2.23) 


0 
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where we deliberately do not explicitly write the integrand in square brackets. Since this 
integral equals 0 for all T E (^“((O, T); C£°- cf. (2.18) - the theorem of du Bois- 
Reymond yields that the integrand should vanish for almost all t £ [0,T] and for //o-almost 
every x. Hence, we obtain 


^ («$,(!)),*!(*)) + j- p («(*«(»)), *!(»)) 


Mdy) 


$*(*) = - VW h {<f> t {x) - $ t (y)) 

n 0 

- V y e (pt($ t (x)), $t(x)) + q[nt](^ t (x), $*(x)). (2.24) 

2.4.3 Step |B] 

The transition to a particle system takes place by substitution of p,ft = Y^i=i d Xifi for po in 
( 2.24[ ). Moreover, in p t and q we replace pt by pf := &t#Po- Note that, after substitution, 
(2.24[) holds pQ-a.e. and should therefore (only) be evaluated at x = Xkp for all k £ {1,..., n}. 


We obtain exactly (2.17). 


2.5 Equations of motion via the route C A]B 
2.5.1 Step [C] 

At the continuum level, deriving the Euler-Lagrange equations resembles considerably what 


was done in Section 2.4.2 Note however that the action as defined in (2.7)--(2.8) is used. In 


(2.7) occurs. The dependence on that is explicitly written down, corresponds 


to the position at which p t is evaluated. However, if is varied, also the function p t itself 
changes. This is somewhat confusing, as this is an implicit, ‘hidden’ dependence of pt on 


the motion mapping <E. However, the exact relation is given by (2.4), which we therefore 


substitute in (2.7). The variational derivative becomes 


S'[§](*) = 



i|i>,(z) + £ *,(x)| 2 - e ( | J(ft ,+ ( ^, )M | .*<W + ^»W 

o d 0 


Po(x) dx dt 


6=0 


cf. (2.20). Some effort is needed to deal with the e-dependence in the Jacobian matrix. We 
refer here to Section 2 of |28| . where the equation of motion is derived from the action, for 
the case where e has no explicit dependence on the spatial coordinate; i.e. e = e(p). The 
determinant of the Jacobian matrix is a polynomial of the entries of that matrix. The basic 
idea in [28] is that the chain rule has to be applied with respect to every element of the 
Jacobian matrix. To avoid having to introduce a considerable amount of extra notation, we 
only state the result of [28( here: 


1 


*.<*) = 


de\ 


$t(x) 


(2.25) 


- ( + pt(&t(x))^(pt($t(x))) ) V#($ t (x)). 
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On the right-hand side the gradient of the pressure P appears, due to the thermodynamic 
relation de/dp = P/p 2 . The reader should note that the notation used in [28] differs sub¬ 
stantially from ours, but that the philosophy of deriving the equations of motion is the same(^] 


If e = e(p,y), and moreover, we include nonconservative forces, then instead of (2.25) we 
obtain 


Og 3^ g \ 

*t(x) = - ( 2-U(T t (x)),<h f (x)) + pi(T t (x))^U(^(.T)),ch t (x))J Vpt(*t{x)) 

- VyG (pt($t(x)), $t(x)) + q[pt]($t(x), $t(x)). (2.26) 


The additional terms follow from similar steps as the ones leading to (2.24). We omit further 
details. Note that, in correspondence with q as introduced before, the dependence on p t in 
square brackets indicates the presence of a nonlocal term; cf. (2.34). In the next step, this 
will become a dependence on the measure pi like before. 

2.5.2 Step [A] 


In this step, we formulate (2.26) in terms of measures. The only place where the mea¬ 
sure pt can be incorporated directly, is in the nonconservative force density. We write 
q[pt](<S>t(x), <3?t(x)) instead of q[p t ](&t(x), All the other occurrences of p t in ( |2.26| ) 

we approximate by pt as defined in (2.9). We obtain 


$t(®) = - ( + /5t($ t (®))|^(/5t(^t(®)),®t(®))) V/5t($t(x)) 

~ V y e (Pt($t(x)), ®t(x)) + q[p t \($ t (x ), $t(x)). (2.27) 


2.5.3 Step [B] 


We evaluate the resulting equation at x = x k p for all k £ {1,..., n} to obtain 


We take pft := Y/a=i m i $ 0 and replace pt by p/ ■= o i R Pt an d Q th a t appear in (2.27 


$t(xk, o) = 

( 3g 3^g \ 

- \ 2-^(pt($t{xk,o)),$t(xk,o)) + Pt{$t(x kfi ))-^(pt{$t(x kfi )),<i> t (x k p))J Vp t ($t(x ki0 )) 

- VyG (p t ($t(x k> o)), $t(xfc,o)) + q[p?]($t(x kt o), $t(x fci0 )), (2.28) 


where each appearance of pt denotes a sum over all particle positions. Namely, 

n 

pt($t(x k> o)) = Y2 m 3 w h($t(x k , 0 ) - ®t(xj, o)), and (2.29) 

3 = 1 
n 

Vp t (® t (x k: 0 )) = ^2 m j o) - ®t(xj, o)). (2.30) 

3 =1 

2 Another interesting observation in [2H] is that the Lagrangian density - when formulated in terms of 
Eulerian coordinates - is just the pressure P. 
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2.6 Comparison of the resulting equations (2.17) and (2.28) 

Procedures 000 and ra yield the same equations of motion, namely ( |2.17| ). As antic¬ 
ipated already in Section 2.2 the equation resulting from Procedure EBB is different; see 
(2.28). This difference between the two final equations arose because we introduced the reg¬ 
ularization via p at different stages. In fact, (2.17) contains an extra regularization in space, 
as we will show now. 


Note that only the term involving de/dp and <9 2 e/<9p 2 is different. In (2.17), we have 

n 

~^2 m i VW h ($ t (xk, o) - 0 )) 


1=1 


dc dc 

— {pt{$t{Xk,o)),®t{Xkfl)) + -Q- (pt(^Ko)),^Ko)) 


while the corresponding part in (2.28) is 

/ d 2 e \ 

- \ 2—(pt($t(xk,o)),®t{xk,o)) + Pt($t(xk,o))-Q^(pt(®t(xk,o)),$t(xk,o)) J Vpi($t(x fej0 )). 

dc 

Note that both of them contain a part — — {pA&tixk o))> ^t(x k o)) ^Pt(®t{x k o))> hence let us 

op ’ ’ 

consider in (2.28) only 

/ p)g d 2 c \ 

- + pt(®t(x k fi))-Q^{pt( , S , t(xk,o)),$t(xkfi))J Vp t (®t(xk, o)) 

= -V fpt(®t(x k ,o))i^(Pt($t(xk,o)),®t(xh,o))') • (2.31) 

To obtain this equality, we have assumed that Vyde/dp = 0; this assumption anticipates the 
choice we make in (2.32). Let us even go back one more step and consider this term before 
the introduction of p, i.e. as in (2.26). To see how this term relates to the corresponding one 
in (2.17), we take the convolution with W k , and proceed as follows: 


- JW h {C - y)Vy (pt(y)^(pt(y),y)] dy = JVyW h {£, - y)pt(y) ^(p t (y),y) dy 

Qt ^t 

f dc 

= - J vw h (£ - $ t (y)) Q^{pt(^t(y)), ®t(y)) Pt(®t{y))\J®t(y)\dy 
n 0 

f dc 

= - J VW4(£ - $t(y)) Q^(pt($t(y)), $t(y)) Po(y) dy. 

Qo 

In the hrst step, we performed integration by parts, with vanishing boundary terms on <912*. 
This is because 12* = suppp* and hence p* vanishes on its boundary. Now replace po(y)dy by 
Po(dy) and approximate p* by p*. Take po : = YH=i m idx i0 and evaluate at £ = $t(xfc,o) and 
obtain 

n „ 

_ 

- } y mj VW h (® t {x k , o) - ®t{xi,o)) w- (p*($*Ko)),$*(x* >0 )). 

i=i ^ 
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This expression exactly appears in (2.17). To summarize: the connection between (2.17) and 


(2.28) is that in the former during the derivation procedure an extra regularization in space 


was introduced for a part of the right-hand side. Note the connection with the following case: 
consider the Frechet derivative £', based on the L 2 inner product, of some energy £ = £(p). 
Define a second energy £ by £{p) := £(Wh * p)- Then £'{p) = Wh * £'(Wh * p), which also 
contains an extra regularization. In this paper we treat a special case of the general energy 
£. 


Note that, if we only consider the part involving de/dp, ( |2. 17 ) is the same as Equation 
(3.8) in [22]. The notation used therein shows the direct dependence on the pressure. In 
Equation (3.5) of [22], the equivalent of (2.28) is given. The reason why (2.17) is the one 


traditionally used in the SPH community is given in [22j : it does conserve linear and angular 
momentum exactly, as opposed to (2.28). Having derived the schemes, we are now able also 
to elaborate on the remark already made in the introduction: (2.28) “requires the gradient of 


the pressure held to be expressed analytically, while the pressure itself does not appear in the 
numerical scheme”. The first part on the right-hand side of (2.28) is - anticipating (2.33) - 

of the form —— — ( p 2 F'(p )) Vp = — — — (P(p)) Vp. Hence we need an analytical expression 


for 


P d P 


p dp 


2.7 Measure-valued formulation 


In Sections |2.3[ 2.4 and 2.5 we derived particle-based schemes. To establish their convergence 
(as n —> oo) we use a measure-valued formulation. Such formulation incorporates both the 


limit and the approximating sequence. Hence, we focus on the measure-formulations (2.24) 


and (2.27), without the specific choice po = Po- Our convergence proof is applicable to a 
class of approximating measures that is much broader than just sums of Dirac deltas. The 


SPH-inspired particle approach is a special case; see Corollary 3.11 


Although (2.24) and (2.27) are different (cf. Section 2.6), we wish to establish the conver¬ 


gence proof for both formulations simultaneously. Hence, we introduce a switching parameter 
6 £ {0,1} to unify both variants in a single equation of motion. First, we assume that e is of 
the form 

e(p,y):=V(y) + F(p), (2.32) 


in agreement with the remark we already made underneath (2.31). Note that de/dp = F' 
and V y e = W. Here, V £ C/(W l : K) describes the portion of potential energy which is due 
to a gravitational or magnetic field and F £ C 2 (M + ;M), where M + := (0,oo) the potential 
energy due to the thermodynamics of the medium under consideration. This decomposition 
of e is typical for an ideal medium, such as a compressible inviscid fluid. Note moreover that 
this is a common modeling assumption in the derivation of the SPH equations for a system 
of particles [22]. We introduce an auxiliary function Fg, 6 £ {0,1}, that is defined by 


1 d 


and F ^p)-=Fip)- 


(2.33) 


We choose q to be of the form 


q[p](y, u) ■■= —rj(y) U+{K * p)(y), 


(2.34) 
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with 77 E C'^(M rf ;M + ) and K E C)J (R d : M d ). The A'-term describes non-local interactions 
within the system, while the 77 -term is a viscous term. We use —77 • u, which is a simplified 
version of the usual viscous term in SPH that (also) involves A Wh * w, see [22]. 


We assign the value 9 = 0 to the formulation in (2.27), and 9 = 1 to (2.24). Both equa¬ 
tions are now simultaneously written as 


(x) = -Fg (pt(Mx))) V/5t($t(x)) - 9 (VW h * [(Fg o p t )p t ])(<S> t (x)) 

- VP ($ t (x)) - 77 ($ t (x)) $ t (x) + (A' * ^)($ t (x)). (2.35) 


Here we use the shorthand notation 


(X7W h * [(Fg o pt)pt]m = / vw h (£ - y)Fo(p t (y)) p t (dy). 


(2.36) 


fit 


In (2.35) we slightly abuse notation, and the equation should be read as follows: whenever 


9 = 0 we disregard the complete term 9 (VIP/, * [(T# 0 Pt)pt])(&t(x)), irrespective of whether 
the convolution term is well-defined, bounded etc. 


Remark 2.1. We emphasize that Ao and F\ are physically different objects in the sense that 
Ao contains all contributions of A to the flow, while F\ only contains part of that influence. 
Hence, although the notation might suggest so, by setting 9 = 1 we are not adding terms. 
We use one function Fg to facilitate the presentation in the sequel. However, Ao and F\ do 
have the same physical dimension and e.g. if A is given by F(p) ~ p K for some k 61\ {0}, 
then both Ao, Ai ~ p K ~ 1 . 

Now we arrive at the central evolution problem we will consider in the rest of this paper. 
Fix a final time T > 0. Let A(M rf ) be the space of probability measures on M d . Assume that 
po E A(M d ) and that there is an tq > 0 such that 


supp/uo C B(r 0 ). (2.37) 

Let vq E Cl (M d ; R d ) and 9 E {0,1} be fixed. We consider the system 

' $t(x) = -Fg (, p t (M x))) V / 5 t ($ t (x)) - 9 ( VW h * [(Fg o p t )p t ])(<S> t (x)) 

-VP ($ t (x)) - 77 ($ t (x)) $ t (x) + (K * Mt)(^t(*)); 

< p t -=W h *p t ; (2.38) 

IH = $t#P 0 ; 

. ^o(x) = X , $o(x) = vo(x), 


for all x E supp/io and all t E (0, A], We remark that this condition implies the one with 
(2.24): that equation is required to hold for almost all t E [0,T] and for /iQ-almost every x. 


Remark 2.2. We might have taken K * (Wh * Pt) for some K, instead of I\ * pt, to comply 
with the pressure term (i.e. the one involving Fg) that only depends on the regidarized density 
pt- We prefer the shorter form K * pt- This choice can be made without loss of generality if 
we take K = K * Wh- 
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Remark 2.3. It is not a priori clear whether the term K * pit is a conservative or a noncon¬ 
servative force density, hence whether it should be part of q or be related to e. Assume there 
is a K such that K(f) = — A^(|£|)£/|£|- Then both ways give the same equations of motion. 
Indeed, if we include the energy density | K * p t in e instead of including K * p t in q , we also 
obtain (2.35). 


3 Convergence 


In this section we introduce some preliminary notions, and summarize the required assump¬ 


tions together with the convergence result (Theorem 3.10). The theorem provides a general 
result, of which the convergence of SPH schemes is a special case; see Corollary 3.11 The 


proof of the theorem in given in Section 3.4 


3.1 Preliminaries 

Fix a constant integer d £ N + . 


Definition 3.1 (Push-forward). The push-forward of a probability measure pi £ V(M. d ) by a 
mapping <I> : —>• W 1 , notation Qjfpi, is defined by 

(4>#/r)(R) :=^($- 1 (i?)) (3.1) 

for all measurable B C M rf . Equivalently, we can define &Jfp as the push-forward of pi by & if 

J fix) {<f>#p)(dx) = J fi<&{x)) n{dx) (3.2) 

R d R d 


for all measurable, bounded functions f on M^. 

Definition 3.2 (Joint representation). A joint representation of two measures fi\, p 2 £ P(M rf ) 
is a measure tt on M. d x M. d such that 

7 t(A x R d ) = p\{A), and 7r(M d x B) = pi 2 {B), (3.3) 

for all A and B in the Borel a-algebra o/M rf . We denote by II(/ui, //s) the set of all joint 
representations of p i and pi 2 ■ Joint representations are also called couplings. 

A useful property of a joint representation 7 r £ II(//i, p 2 ) is that for each i = 1,2 


/ 


f(xi)Tr(dx 1 ,dx 2 )= f(x)pi(dx) 


(3.4) 


for all measurable, bounded functions / on M rf . In fact, this is an alternative definition. 


Definition 3.3 (Wasserstein distance). The Wasserstein distance between two probability 
measures pi,p 2 £ P(M d ) is defined as 


1, M2) 


inf 


x — y | ir(dx, dy). 


(3.5) 
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Note that, to be more precise, we should call this the 1 -Wasserstein distance , as a special 
case of the p-Wasserstein distance for which the cost function \x — y\ p is used in the integral. 
The 1-Wasserstein distance is usually written as W\ , but we will stick to W to avoid confusion 
with the smoothing function W/The particular choice p = 1 is made because it is compatible 
with the Lipschitz properties of the functions and the motion mapping that we use. This is 
what Section [3] hinges on. For an exposition on the Wasserstein distance and the related 
concept of optimal transport , we refer to [29] and EDI- 

3.2 Assumptions 

Throughout the paper, we assume the following: 

Assumption 3.4. The functions V, q and K satisfy V G r? G C^(M rf ;M + ) and 

K G cf(R d -,R d ). 

Remark 3.5. Note in particular that the above assumption implies that VU and K are 
Lipschitz continuous. We denote their Lipschitz constants by | W|l and \K\l, respectively. 

For Fg and Wh we have requirements that depend on the value of 9. Recall that 

M + := ( 0 , oo) 

and define 

Mp" := [ 0 , oo). 


Assumption 3.6. The function Wh G is even and satisfies f Rd Wh(x) dx = 1. 

Assumption 3.7 (The case 9 = 0). We require that Fq G Moreover, we assume 

that there is a constant M\ > 0 such that for all p, G R(M d ) 

sup \F 0 ((W h * p)(x))X7(W h * p)(x)\ ^ M 1 . (3.6) 

ccEsupp fii 

If 9 = 0, we define M 2 , M 3 > 0 such that 

sup l-Fo('u)! ^ M 2 , and (3-7) 

uGUr,w h 

sup |F»| <M 3 , (3.8) 

ueUr,w h 

where 


Ut,w h ■= G M + : f ^inf^^ W h ) IjlFhlloo | , and (3.9) 

r(T) := r 0 + T||u 0 ||oo + \t 2 (||W|U + M l + ||(3.10) 


cf. (3.25). Under Assumption 3.7 Fq may have singularities at the origin, but only if Wh is 


strictly positive everywhere in B(2r(T)). Such Fq and Wh are used in [7]; see also Section 3.5 


If 9 = 1 we need the following assumption: 
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Assumption 3.8 (The case 9 = 1). We assume that F\ E C' 1 (Mq";M). 

For 9 = 1, let M 2 , M 3 > 0 be such that 

sup |Fi(tt)| ^ M 2 , and (3-11) 

ue[o, WW^} 

sup |F((ri)| ^ M 3 . (3-12) 

ue[o,||W/i||oo] 

and define M\ := 2M 2 HVW/Jloo- 

In both cases 9 = 0 and 6 = 1, we use the same letters for the constants, to ease nota¬ 
tion in the sequel. 


Remark 3.9. The upper bound in (3.6) is needed to get an a priori bound on the propagation 


speed in Lemma 3.12. Consequently, we can restrict ourselves to measures with bounded 


support afterwards; cf. Corollary 3.13 To achieve Lemma 3.12 if 6 = 1, we need Assumption 


3.8, which does not allow for singularities in F\ around zero. 


We demonstrate now why a weaker assumption for F\, resembling (3.6) is not feasible. Assume 


that F\ (p) := p a with a E (—1,0). This is the case also considered in [?]. To bound the first 
term on the right-hand side of (2.35), in [7] it is assumed that |VIF/j(£)| ^ c\Wh{0\~ a for 
some c > 0. We would need an estimate on 


sup | (yW h * [(Fi o (W h * At)) 

ccEsupp fj, 


m])0)| 


(3.13) 


Let Wh be strictly positive everywhere. Since Wh G L 1 (M d ), lim^oo Wh{Q = 0. Let Wh 
satisfy the aforementioned condition |VtFft,(£)| ^ c\Wh{C)\~ a - Then also lim^oo |VW/ l (^)| = 
0. Under these (not very strict) conditions one can show that (3.13) is unbounded; to see this, 


use e.g. the sequence of measures {p K ) K ^+ defined by p K := (<5_ ree] + 5^ + <5( K+1 ) ei )/3, where 
e L is the first unit vector in M d . Note in particular that (3. 13|) is unbounded for a Gaussian 


Wh- The Gaussian however is one of the standard choices for Wh that we do want to allow 
for. 


3.3 Main convergence result 

Let {a*o }ngN C - P(M‘ i ), and assume that 

suppAto C B( r o) for all n E N, 


(3.14) 


where ro > 0 is the same constant as in (2.37). For each n E N we associate to the measure 
Pq a system of equations analogous to (2.38[): 


$?(*) = ~Fe (/5?(^(s))) Vp?($?(s)) - 9 (V^ * [(F e o ft) 

-VU ($?(*)) - 7/($?(®)) $?(s) + (K * 

ft := Wh * ft; (3.15) 

(ft = 

, ^0 ( x ) = X 1 $&(*) = ^0(3;), 
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for all x E supp/Z q and all t E [0, T]. Note that the only difference with (2.38) lies in the 
initial distribution gfi versus /zq; the initial velocity vq is the same. 


For any r > 0, define V r (M . d ) := {/z E V(R d ) : supp n C B(r)}. We also define A as the 
space of all functions from supp po to C 2 ([0, T];R d ). 

The main result of the present paper is the following. 


Theorem 3.10. Assume that vo E 
moreover (depending on the value 
sequence {/Zq} C TV 0 (M d ) such that 


cl( 


pd. Tn>d\ 


and that Assumptions 3-4 and 3.6 hold. Let 
of 8) Assumption 3.1 or 3.8 he satisfied, and take the 


W(/z^ 0 )™0, 


(3.16) 


for some po E P ro (M“). Then: 


1. there is a unique pair (/z, <F) E C([0, T]; V r r r )(W i )) x A that satisfies (2.38); 


2. if, for all n E N, the pair (/z n ,<h n ) E C([0,T];x A is a solution of (3.15), 
then 

sup W(^)™ 0 . (3.17) 

te[o,T] 


As a corollary, we obtain the following convergence of the SPH scheme with n particles. 


Corollary 3.11. Fix 6 E {0,1}. For each n E N + , let p]f := Jf 7 j=i m j^j o e TV 0 (R d ) /or 
some {mj} 7 j =1 C M + such that = h anc ^ / or some {xj t o}”=i C B(ro). Assume that 

W(/Zo,/zo) ”2128 o for some /zo E T-V^R**). T/ien f/ie discrete measure pf = Ylk=i m k^ t (x k 0 ) 
associated to the particle scheme defined for each k E {1, ..., n} fry: 


$t(zfc,o) = ~ v ^(^i( x fc,o) - ^z(^i,o)) [-fe (Pt($t(xk, o))) + 0F e (pt($t(xi,o)))] 


%— 1 


- W ($i(x fc ,o)) - J?($t(aofc,o)) $z(xfc,o) + (K * Pi)($ t (xk, o)), (3.18) 


converges to the solution pt of (2.38) in the following sense: 


sup w(i$,ih) 0 . 

te[o,T] 


(3.19) 


3.4 Proof of the main convergence theorem 


Before proving the main result, Theorem 3.10, we need two auxiliary lemmas concerning the 
properties of the motion mapping The first lemma is an upper estimate for 
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Lemma 3.12. Let Assumptions \3.4\ 3.6 and \ 3. 7| or 3.8 (depending on the value of 9) he 


satisfied. Then for any given g G C([0,T]; V(R d )) the mapping in (2.35), completed with 
pt := Wh * pt, 3>oOe) = x and 3>o(x) = vq{x), satisfies 


|$i(x)| < |x| + 
for all x G supp and all t G [0, T]. 


oo + \t 2 (Mx + IIvyIloo + Halloo), 


(3.20) 


Proof. For p fixed, and for each x G supp^o, the ODE (2.35) is well-posed on [0,T], given 
the assumptions on V, ig Fg, Wh and K, and the fact that p is continuous in time. The 
well-posedness follows from the Picard-Lindelof Theorem; further details are omitted. 

Using an integrating factor H(t) := exp ^ JqP{^t{x)) dr'j, we deduce from (2.35) that 


|<f>t(x)| ^|<FoOe)| + \vo(x)t\ + 


t S 

r i 


J H(s)J 
0 0 


t S 


^\x\ + t 


+ 


0 0 


H(r) 
H(s ) 


H{r) (<£ r (x) + r/{& r {x)) <h r (x)) dr ds 


dr ds. 


<f> r (x) + rj($ r (x)) 6 r (x) 


Since i] is a positive function and hence 0 ^ H{r)/H{s ) ^ 1 in the inner integral, it follows 
that 


t S 


|«F t (x)| ^ \x\ + t 11 Vq11 oo + 


-X7V(Mx)) + (1C * Vr)(Mx)) 


0 0 


- Fg {{W h * Pr)(Mx))) V(VF h * Pr)(Mx)) 

- 9 (VWh * [{Fg o {Wh * p r ))pr]){Mx))) 


drds. (3.21) 


In the case 0 = 0, the following estimate holds due to Assumption 3.7 


| Fg {{W h * Pr)(Mx))) V{W h * g r )(<S>r(x)) + 0 (VWh * [{Fg o {Wh * Pr))p r ]){Mx))) \ 

= |F 0 {{W h * PrfiMx))) V{W h * Pr){Mx))\ ^ M x . (3.22) 

Note that for any p G P(M d ) it holds that \\Wh * /i||oo ^ 11 UUft. 11 oo ■ Hence, in the case 0 = 1: 

| Fg {{Wh * Pr)(Mx))) V(W h * Pr){Mx)) + & {VW h * [{Fg o {W h * p r ))p r ]){<S> r {x)))\ 

^ M 2 HVWfclloo + IIVWfclloo M 2 = M\, (3.23) 


where the bounds of Assumption 3.8 are used. 


A combination of (3.21), (3.22) and (3.23) yields that for each 0 G {0,1} 

t S 


|$i(x)| ^ |x| +t Hwolloo + J J (|| VU||oo + Halloo + Mi) drds 

o o 

holds for all x G supp and t G [0, T ], from which the statement of the lemma follows. □ 
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ro \ 


3.6 and 


Corollary 3.13. Let fiQ E V- 

on the value of 9) be satisfied. Then any solution of (2.38) must satisfy 


and let Assumptions 3-4 


3.1 or 


3.8 (depending 


supp fit C B(r(t)), 


(3.24) 


for each t E [0, T], where 

r(t) := ro + i||vo||oo + ^ i? (Mi + ||VV 11oo + || K 11oo)• (3.25) 

The next lemma provides a Lipschitz-like estimate on 4?f. 

Lemma 3.14. Let v l ,v 2 E C([0,Tj ;be given. Consider the motion mappings 
corresponding to v l (i = 1,2): 

*i*(0 =-Fo ((W h * ^)(<(0)) V(W), * (3.26) 

- 9 (vw h * rn o (W h * 4)>i])(<^)) 

- W«(£)) - + (K * 

for all £ E supping an d t E [0,T], with initial conditions 4>q 1 (£) = <1>q ! (£) = uo(£). Then, 

for all t E [0, T], x E suppz/Q and y E suppr^, it holds that 

Wt ( x ) “ $f(y )I < (! + t Halloo) \x-y\+t |u 0 (x) - vo(y)\+ 

t t 

+ J[Mi(t-s) + \\ri\\ 00 \\§f(x)-$>f(y)\ds + M 5 J(t - s)W(v), v 2 s ) ds, (3.27) 

0 0 

where 


M 4 :=\VV\ L + (l + 9)M 2 \\D 2 W h \\ 00 +M 3 \\VW h \\ 2 00 + \K\ L , and (3.28) 

M 5 := (1 + 9)M 2 ||H 2 lT /t || 0O + (1 + 9)M 3 ||VW fc |& + \K\ L . (3.29) 


Proof. Note that, by the Fubini’s theorem, for any integrable function /, we have 


t r t t t 

f(s)dsdr = J J f(s)drds = J (t — s)f(s) ds. 

0 0 0 s o 


(3.30) 
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Integration of (3.26) in time together with (3.30) yields that 


M (*) - W (y)l <1* -y\ + t M®) - v 0 { y )\ 


(3.31) 


Z 

+ J(t- s)\VV(^ s \x)) - W(<hf (y))| ds 


+ 


+ 


t r 


v(*s (x)Ws (x)-T m{y)W a (y) ds dr 


0 0 
t 


(t - S ) | Fg (iW h * ^)(<&f (*))) V{W h * ul)(^(x)) 
-F e [{W h * v 2 s ){*f{y))) V{W h * v 2 a ){*?(y)) 


ds 


+ ej(t- 8 ) (yw h * [{F e o {w h * ^)K])(<^ (*)) 

o 

-{VW h * [{F e o {W h * v 2 s )y s mf{y)) 


ds 


+ J{t~s) {K*^m {x))-{K*vj){& a {y)) 
o 


ds. 


Furthermore, we have 


Z Z 

[ {t - S )|Vy(^(x)) - W(<hf (y))| ds < |VV|l J{t- S )|^ 1 (x) - (y )\ ds, (3.32) 


and 


t r 


t r 

-J^V) 


[ [ 0))$f (*) - (y))$f (y) ds dr 

= 

HI 

/ rj{z) dz 

ds dr 

J J 

0 0 


0 0 

J 

®f(y) 



t 

p 

y ( x ) x 

r r 


I 

/ rj{z) dz — rj{z) dz 

dr 

0 

ff(y) y 



t 

H\ 00 J\^ 1 {x)-^f{y)\dr+H\ 00 t\x 

o 


— 2 / 1 - 


(3.33) 


Regarding the term involving Fg on the third and fourth line of (3.31), we proceed as follows 


|Fg {{W h * ^)(£i)) V(W h * ^)(£i) - Fg {{W h * v 2 ){&)) V(W h * v 2 ){&)\ (3.34) 

< |Fe {{W h * ^)(6))| IV(w fc * ^)(6) - v {W h * ^)(6)| 

+1 Fg {{W h * ^)(6))| I v{w h * ^)(6) - V{W h * ^ 2 )(6)| 

+ I Fg {{W h * ^)(<£i)) - Fe {{Wh * ^)(6))| I V(w ft * ^ 2 )(6)| 

+ j Fg {{W h * 1/^(6)) - Fe {{W h * I' 2 )(6))| |V(Wfc * *?)(&) j • 
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We only consider 6 E supp and 6 G suppzz 2 . This implies that 6)6 E S(r(T)). For each 
i, j E {1)2}, we have the following estimates: 

(Wfc *!/*)(&) < Halloo, 

(since z/* is a probability measure), and 

(W h * z4)(6) > inf W h (6 - z) = inf W h . 

K a)K ™' 5i) *eB(r(T)) h ^ J ’ B(2r(T)) 


Thus we get (Wh * ) (6) G Ux,w h - We proceed with the estimation of ( |3.34 ): 

|*b ((w A * ^)(6)) V(w h * ^)(£i) - Fg (( w h * z/ 2 )(6)) v(W), * ^ 2 )(6)| 

< M 2 IIZD 2 VF/j.||oo |6 - 61 + m 2 I V(W h * ^)(6) - V(W h * i/ s 2 )(6)| 

+ M 3 ||VW h ||L I£i -&| + Af 3 ||VWft|| 0O |(W /l ^ s 1 )(6)-(W /l ^ 2 )(6)|, (3.35) 

where we used that ||V(W fe * zz 2 )||oo < ||VW^||oo, \\D 2 (W h * Vg)\\oo ^ \\D 2 W h \\ 00 and the fact 
that \tl >\l = HV^Iloo for any differentiable function ip. Note that: 


\(W h * i/g)(&) - (Wh * zz 2 )(6)| = 


Wft (6 - 2 ) ^(<fe) - Jw h (& - w) V 2 s (dw) 


J(W h (& - z) - W h (& ~ w)) Tc s (dz, dw) 


< J\w h (& ~ z) - Wft (6 - w)| 7 r s (dz, dw) 
^ IIVW/1II00 J\z-w\ TT s (dz, dw), 


where tt s G II(zz},z/ 2 ) is arbitrary and the second equality follows from (3.4). By minimizing 
over all couplings in II(z/}, u 2 ), we obtain 


\(W h *V l M2 )-(W h *V 2 M2)\ SS HVWftllooW^ 1 ,^ 2 ). 


(3.36) 


We stress that the bound (3.36) is independent of the choice of 6>6- Analogously, we have 
\V(W h * ^)(6) - V(W h * ^ 2 )(6)| P'WfcllooW^ 1 ,^). (3.37) 

It follows that 

| Fg ((Wh * ^)(6)) v(w h * ^)(6) - F e ((Wh * ^)(6)) V(w h * *?)(&) | 

< (M 2 \\D 2 W h \\ 00 + M 3 \\VW h \\ 2 00 ) (|6~6l+W(^,zz 2 )). (3.38) 


If 9 = 1, similar estimates as in the first term on the right-hand side of ( 3.35[ ), and as in 
( |3.36[ ) and (3.37) yield 


\(VWh * [(Fg o (Wh * ^))^])(6) - (Vlbfe * ((Fg o (Wh * ^))^])(6)| 

< M 2 \\D 2 W h \\ 00 \ii-i 2 \ + (M 2 \\D 2 W h \\ 00 + M z \\VW h \\l 0 ) W(i/ a \i/ s 2 ). (3.39) 
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(3.40) 


The last term in (3.31) we treat as follows: 

I (k * vim) - (k * u 2 s m) | <\{k* vim) - & * vim) I 

+ \(K*vim)~(K*v 2 s m)\ 
da- 61 +w(vi,m). 

The estimate of the second term is obtained like in ( |3.36 ). 

We combine flOI] ), ( [3732] ), ( [3733] ), ( [3738] ), ( fT39] ) and ( f330| ) to get 
\$f (®) - ®f(v)\ < (! + * Ihlloo) \x-y\ + t \v 0 {x) - vo(y)\ 

t t 

+ J[M 4 (t-s) + \\ri\\ 00 ]\$>f(x)-<i>f(y)\ds + M 5 J(t. - s)W(vl,v 2 s ) ds, 

o o 

with M 4 and M 5 as defined in the statement of the lemma. □ 


Now we have all ingredients to prove the main theorem, Theorem 3.10 

Proof of Part [Q of Theorem 13.101 

If M 5 = 0 then the well-posedness of (2.38) is straightforward. In that case, Fgo{Wh*yt) = 0 
on supp m for all t and moreover K must be constant, so the first equation in (2.38) is inde¬ 
pendent of yif The Picard-Lindelof Theorem guarantees, for each x E supp/ro> existence and 
uniqueness of the motion mapping (as mentioned before). The solution (yt)o^.t^T is uniquely 
defined by the push-forward yt = 


If M 5 0, the well-posedness proof is based on a fixed-point argument (Banach’s Fixed 
Point Theorem). Let T > 0 be fixed. Choose N E N + large enough, such that T* := T/N 
satisfies 

k t * := i (T *) 2 M 5 exp Urjl^ T* + i M 4 (T*) 2 ^ < 1. (3.41) 

Let j E {l,...,iV} be hxed. Suppose that yip E P(M d ) and vP E C^(M d ;M d ) are given. 
Consider a mapping : v e -> y := J-^\v) from 


Cj := {v E C([0,T*];V r{j mR d )) : Ht=o = (3-42) 

to itself, defined by 

IH = for all t E [0, T*], (3.43) 

where the motion mapping v : supp yip —> C 2 ([0, T*]; M d ) is the solution to the following 
ODE 




[*F } n*) = Vpt([$W)Ux)) 

-0 (VW h * [(Fg o pt)vt])([*?nx)) 

-VP ([^V(x)) - T] ([^H®)) [$Pnx) + (K * ^)([$?Y(x)); 

Pt '■= W h * vy 

%(x) = x, 4>£(x) = vp\x). 

(3.44) 
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The space Cj is complete for arbitrary j E {1,..., N} due to Theorem A.4 in Appendix [Aj 
Note that a fixed point /jP^ of this mapping together with the corresponding motion mapping 
is a solution of (2.38) on [0, T*] with initial data /ijp and Vq. We create a hierarchy 


of the mappings for j = 1,..., N by defining 


and w 0 := vq. Such definition only makes sense if mapping j actually has a unique fixed 
point and thus ppi and T)-* are well-defined. Moreover, we are aware of the fact that we 
have only defined 'Cq' 7+1 ^ on the support of This is however sufficient. If we in¬ 

sist, we might just define it to be zero outside. In view of the to be constructed hierarchy, 
supp/x^ C B(r((j — 1 )T*)) should be satisfied for each j. 


A*o 


0+1) ._„()) „(!)._ 


— firpt, /ig flQ, Vq 


L + 1 ) <jjO') 


For any v E Cj the image /i = J-- 3 ^(u) exists, and actually is an element of Cj. Well-posedness 
of the motion mapping (for given u and for each x E supp/Xg^) follows from Picard-Lindelof 
(see before) and guarantees the existence and uniqueness of /x. 

The support of the image measure, supp/x?, is contained in a ball of radius 


r(jT*) = r 0 +jT* 


+ ±(jT*) 2 (M 1 + || W||c 


+ \\K\ 


(3.45) 


This is easily checked by use of (|3.20 ) and a recursive relation involving ||[<F 


,()+!) I 

iTo I 
ment of C 


k(?)i 


for each j E {!,..., N — 1}. Thus, the image /x of our mapping is an ele- 


3- 


Consider two measures u 1 ,^ 2 
p 2 := F^\v 2 ). Let 7To E II (/x^ 
fine 7 T t E by 


E Cj and their corresponding images /x 1 := and 


, /Xq ) be arbitrary. For an arbitrary fixed t E [0, T*] de- 


n t ■■= [$ 




>[* 


U) i 


#^ 0 - 


(3.46) 


Note that this ir t is indeed a joint representation of /xj and /j 2 for each t. We drop the 
dependence on j of 7To, A^. 1 , yu 2 and n t since no ambiguity appears. By definition of the push- 


forward and of the Wasserstein distance (see Definitions 3.1 and 3.3), we have 


I \z-w\ir t (dz,dw) = I [4>| j) ] 


'\x) ~ vr 0 (dx,dy) 


(3.47) 


holds for each t E [0,T ]. Applied to (3.27), a version of Gronwall s Lemma yields that for 

(?) - 

each 1 , 1/6 supp /ig 

[#r 1 (x)-[®lV( S ) 




(1 + t IMloo) \x-y\+t |^ j) (x) - uJ j) (y)| 


+ M 5 /(f - s)Wpl, v 2 ) ds 


exp 


|oo t + - M4 C 


(3.48) 


We remark that Gronwall’s Lemma may be applied because the term |ug^(x) — Ug^(y)| is 
bounded and s e -> (t — s)Wpl, u 2 ) is bounded and continuous. The former can be shown by 
using estimates similar to those in the proof of Lemma 3.12 The boundedness of the latter is 
trivial, 
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(t — s)yV(ul,u'g) < T* sup se r 0 T »] , Vg); while the continuity of s ^ W{v].,vj) follows 

from the triangle inequality. Indeed, since 


m^ s y s ) - i < 

it implies that lim s ^ So W(z^, v 2 s ) = W{v] , v 2 s ) if W^J, ^ ) ->• 0 and W(z^, ^ ) -> 0. 


Now we combine (3.47) and (3.48), and obtain 




(1 + t Halloo) J\x - y\n 0 (dx,dy) + t J\v^\x) - v^\y)\Tr 0 (dx,dy) 


+ M 5 J(t- s)W{vl,Vg)ds 

o 


exp 


\oot+-M A t ). (3.49) 


The integral with respect to Tro(dx,dy) disappeared for the third term inside the square 
brackets, since this term is independent of x and y, and moreover f no(dx,dy) = 1. Now we 
take 

7T 0 := (I® I)#Ho\ 

(j) 

which is the measure concentrated on the diagonal x = y with marginals both jUq . With 
some abuse of notation it can also be written as 

ir 0 (dx,dy) := 8{x - y)$\dy). 

For this choice of 7To, we have that 


j\x — y\ vro(dx, dy) = 0, and 


\vo\x ) - v^\y) I ir 0 (dx,dy) = 0. 


Therefore, only the third term in square brackets on the right-hand side of (3.49) remains. 


Since 


t 

J(t- s)W(vl,v 2 )ds < 
0 


t 


sup W(z/i, u 2 ) / (t — s) ds 

sG[0,t] J 


\t 2 sup W(vl,v 2 ), 
z se[o,t] 


we obtain 

W(/il, l-if) \t 2 M 5 exp (Halloo t + \ M 4 tA sup W(z/', v 2 s ). 

1 \ z J se[o,t] 

Finally, we take the supremum over t 6 [0, T*\. 

sup (T*) 2 Ms exp f H^loo T* + )- M 4 (T*) 2 ] sup W(^\ vj). 

te[o,T*] z \ z / te[o,T*] 
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By the specific choice of T*, J 7< ' J> is a contraction mapping for each j, since 


sup W (/xj, ) < k t * sup W(vl v ?), 

te[0,T*] t£[0,T*] 


where kt* < 1 by assumption; cf. (3.41). As mentioned before, the space Cj is complete for 
each j due to Theorem A.4 in Appendix [Aj Banach’s Fixed Point Theorem then guarantees 
the existence of a unique fixed point of for each j. 

Having the construction of (//■?), $(■?)) f or j = \ : ... ^ N, we define a couple (/x, 4>) of a measure 
and a motion mapping as follows 


*t) := (/4-(j-1)T* ’ t G (O' - l)r*> jT*], 


di) 


(3.50) 


for j G {1,... ,iV}. 

By our construction (/x, 3>) G C([0,T]-,V r (T) 
initial data /xq and vq. 


x A and it uniquely satisfies (2.38) with 


Proof of Part [2] of Theorem 13.101 

Note that Part [I] implies that for each initial measure /xo and (for each n G N) there is a 
corresponding unique solution (/x, <f>), <h"), respectively. Fix n G N and let 7To G n(/XQ , /xo) 

be arbitrary. We use ( |3.27[ ), taking u 1 = y n and u 2 = /x. Thus <h !/i = <f> n and 4W = <f>. First 
of all, we estimate 


|u 0 (x) - v 0 (y )| ^ ||Vu 0 1|oo |a; - y|, 


(3.51) 


for all x G supp /Xq and all y G supp /no • This is possibl^J 
defined on the whole of and has bounded derivative, 
integrating (3.27) against TTo(dx,dy), we obtain 


since vq G C£(M d ; M d ) is given, is 
Using this Lipschitz estimate and 


J mx) - $t(y)| n 0 (dx, dy) < (1 +t (|| Vv 0 ||oc + IMU)) J \x - y\ ir 0 (dx, dy) 

t t 

+ J [M 4 (t — s) + llr/Hoo] J\®s(x) - <& s (y)\TT 0 (dx,dy)ds + M 5 J{t - s)W(/x", /x s ) ds, 

0 0 

(3.52) 

where we used that the last term is independent of x and y , and the fact that 7To is a probability 
measure on xlf If we define tt s G n(/x”, /x s ) as 

7f s :=($?, <& s )#tt 0 (3.53) 

for each s G [0, T], then we have, analogously to ( |3.47[ ), the following: 

W(/x™,/x s )^ J\z - w| n s (dz,dw) = J\^{x) - $ s (y)| ir 0 (dx,dy). (3.54) 


3 Note that this estimate was not possible in ( |3.48| ), since Vq' 1 is part of the solution and only defined on 
supp/io^- In general, Vvg^ might not even be defined. 
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We substitute this estimate for W(Ms, y s ) in the right-hand side of (3.52) and apply Gronwall’s 
Lemma to obtain 


\$t(x) ~ $t{y)\ 7 T 0 (dx,dy) 




(1 + t (II VuqIIoo + IMloo)) \x-y\ 7T 0 (dx, dy) 


Together, (3.54) and (3.55) yield 


exp 


|oo t + — (-M4 + M5) t 2 


(3.55) 


VW>A*t) 




(1 + t (||Vu 0 1 |00 T 11 77 11oo)) / \x - y\n 0 (dx,dy) 


exp 


100 t + — (M4 + M5) f 2 


We take the infimum over ttq E II(/Xq,^o) on the right-hand side: 


>V(Mi\Mt) ^ (1 + * (l|Vu 0 ||oo + Ik/Hoo)) exp 


,t + - (M 4 + M 5 ) t 2 ) W(/$, mo). (3.56) 


Finally, we take the supremum over t E [0, T] on both sides of the inequality and obtain 


sup W(nt,nt) ^ (1 + T (||Vu 0 1 |00 + IImIIoo)) exp 

ie[0,T] 


,T+-(M 4 + M 5 )T z ) w(m^,mo). 


Hence yV(Mo,Mo) n —? 0 implies sup tg [ 0 T ] W(nf,y,t) n —> 0. This finishes the proof. 


3.5 Discussion on Assumptions 3.7 and 3.8, and the condition (3.16) 


We comment on the assumptions needed for the main theorem, Theorem 3.10 


Assumptions on Fg and W^: We remark here that in [7] only 8 = 0 is used, and further¬ 
more VH = 0, 77 = 0 and K = 0. All possible Fq and Wh treated in [7| satisfy Assumption 


1. Fq(u) = u a , for a ^ 0, satisfies the assumptions for all choices of Wh E C 2 (M d ;] 


u) h 


2 . Fq(u) = u a , for — 1 < a < 0, satisfies the assumptions if Wh is an element of C%(W d -, R + ) 
and satisfies the extra condition |VH4,(x)| ^ c\Wh(x)\~ a for all x, for some constant 
c > 0 . 


We remark that the class of admissible pairs (Fq, Wh) covered by Assumption 3.7 


is more 


general than in [7], where only To of the form Fq(u) = u a is treated. For instance, in our 
work any Fq E C^(M + ; M + ) (bounded and with bounded derivative) is allowed in combination 
with an arbitrary Wh E C 2 ( 




Assumption (3.16) on convergence of initial data: Given the initial probability measure mo 
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supported in the ball B(tq), we demonstrate here two ways of constructing an approximating 
sequence of measures (Mo)neN+- 


The first way of constructing /Xq is deterministic and has been used in j3j. For simplicity 
of presentation, we assume d = 1 and suppMo C [0,1]. For each n E N + , define 

n 

juq := rrii6±_ j_, (3.57) 

* ^ n 2 n 

i —1 


where rrii : = 


fio(dx), for each i = 1,... ,n — 1, and m n := 


Mo (dx). 


It follows that Yli m i = I Mo (dx) = 1 and fi q E 'P(M). Dehne a map T : [0,1] —> ^ : 

1 < i < n} by \k(x) ^ if < x < - and 'F(l) := 1 — For every measurable and 

bounded function /, defined on [0,1] it holds that 


f f( x )noi dx ) 

[ 0 . 1 ] 


i =1 x 7 

n-1 /. / ■ 

2 J dx ) / 



+ J no(dx)f(l 


1 

2 n 


= f f{'&{x))fi‘o(dx). 

[ 0 , 1 ] 


Hence, Mo = ^^Mo- Note that |z — 'F(x)| ^ ^ for every x E [0,1]. Therefore, 


W(mo,Mo)< f \x--$(x)\no(dx) < ^ j no(dx) = ^, 


[o,i] 


[o,i] 


where we obtain the first inequality by taking 7r E n(Mo,Mo) to be tt := (J <8 > ^)#Mo- This 
implies that W(mo,Mo) n —? 0. 

This procedure generalizes to the case d > 1 (but with more involved notation). Let 
supp /ro C [0, l] d and let n E : & E N + }. Dividing the hypercube [0, l] rf into n equal 
subcubes, we obtain analogously that the convergence rate is 0(1/ y/n). 


The second way of constructing /Xq is probabilistic and is based on the law of large num¬ 
bers as already pointed out in (7]. Suppose that the points 7Q, i = 1,... ,n are independent 
identically distributed random variables with the same distribution /xq E V(B(ro)). Let /Xq 
be the empirical measure, defined by 


1 11 

^I>‘- 

1=1 

Note that in fact there is an underlying probability space D and Xi : D —> B(ro). Hence Mo 
is, strictly speaking, not a mere probability measure, but a mapping from D to V(B(ro))\ 
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i.e. /Iq : VI —> V(B(ro)). According to [TOJ, Theorem 11.4.1, the sequence (//q) converges 
almost surely to /jLq. This implies that for almost every realization x±,X 2 , ■ ■ ■ the corresponding 
sequence of measures (/2q) C V(B(ro)) given by /2 ft := 1 converges in the narrow 

topology to n 0 : 

j f(x)no(dx ) ->• J f(x)no(dx), for all / G C b (B(r 0 )). 

B(r 0 ) B(r 0 ) 


The term ‘almost every realization’ refers to the fact that the set (in fi) on which the narrow 
convergence does not hold, has zero probability (with respect to the probability distribution 
on V2). In layman’s terms, this means that if we draw a random sample X\,X 2 , ■ ■., it is ‘un¬ 
likely’ that the corresponding sequence (/2 q) does not converge narrowly. 

Assume that our random sample did yield such narrowly converging sequence (JIq). Since 
all £Lq are probability measures on a bounded domain B(ro), their first moments are uni¬ 
formly integrable (i.e. uniformly in n). Thus, Theorem 7.1.5 in [lj implies that 

o) n -^0. (3.58) 


3.6 Numerical illustration 


We illustrate the theoretical convergence result of Theorem 3.10| by two numerical examples. 
The first one involves only the hydrodynamical force, as described by the first term on the 



equation, via a leapfrog algorithm with a constant time step. The leapfrog algorithm is a 
second-order symplectic integrator with the property of preserving the momentum of the 
system. The Gaussian function, defined by 


W h (x) : = 


1 


-\x\yh 2 


kyjTT 


(3.59) 


for all x € M, is used for the regularization of the mass measure in the one-dimensional case. 
For d = 2, the cubic Wendland function is used, whence for all x G M 2 : 


W h (x) : = 


'|(1 + 2,\x\/2h){2 - \x\/h) 3 , \x\<2h, 

0, \x\ > 2 h. 


(3.60) 


These choices are made to illustrate that we can handle both bounded and unbounded support 
of W h . 

In order for the regularized equations of hydrodynamics to approximate the real physics well, 
h should be sufficiently small. Let Vo denote a representative volume assigned to each particle 
based on the initial configuration. In a bounded domain, typically Vo scales as Vo ~ 1/n. It 
is common practice to achieve "h sufficiently small” by taking h = e -^Vo, with parameter 
1.2 ^ e ^ 1.5, cf. |22j . However, the convergence result in Theorem 3.10 holds for h fixed, 
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cl +1 



7 8 9 


Figure 3.1: For 7 = 1 , density p at final time T = 1, with n = 2 9 particles and convergence 
C 1 . Red and blue plots refer to the schemes for 6 = 0 and 9 = 1, respectively. Upper plots 
present results for h = 1 fixed, while lower plots to variable h = 1.5 Vq. 
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Figure 3.2: For 7 = 2 , density p at final time T = 1, with re = 2 9 particles and convergence 
C 1 . Red and blue plots refer to the schemes for 6 = 0 and 9 = 1, respectively. Upper plots 
present results for h = 1 fixed, while lower plots to variable h = 1.5 Vq. 


and the dependence of h on n is not investigated. Numerically, we investigate both cases. 
That is, we take both h = 1 fixed and h = 1.5 \/Vq, which hence varies with the number 
of particles. We assume that the initial measure po has a density po such that po(x ) = 1 
for all x € [0, l] d and po(x) = 0 otherwise. We construct the measure Pq, corresponding to 
the re-particle approximation, according to (3.57) or its d-dimensional counterpart. Hence, 
the initial particle configuration is realized for d = 1 by equipartitioning the initial domain 
[0,1] into re volumes. For the two-dimensional examples, the initial domain is the square 
[0, l ] 2 C M 2 and particles are placed in the center of each square incremental volume partition 
Vo- Masses are assigned as rep = po(xi ) Vo for each i = 1,..., re. Note that in case of more 
complicated initial domains, an equipartitioning may be obtained with a centroidal Voronoi 
tessellation. 

As argued underneath ( 3.57[ ), the sequence (^o) ne N+ constructed in this way converges to po 
and the convergence rate is 0(1/\fn). Hence, the corresponding solutions (p n ) n eN+ converge 
at the same rate; see the last lines of the proof of Part [2] of Theorem 3.10[ 


The hydrodynamical problem considers the spontaneous expansion of a gas cloud until time 
T = 1 , governed by the equation of state P(p) = /Cp 7 , where /C = 1 is a parameter and 
7 the so-called polytropic exponent. We recall that P relates to e via de/dp = P/ p 2 . In 
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dimension d = 1, we examine the cases 7 E {1, 2, 7}, using a constant time step At = 10 3 , 
and the Gaussian function. Note that the case 7 = 1 is not covered by the convergence proof 


(cf. Assumption 3.7 and Section 3.5). It is a limit case (the proof does hold for any 7 > 1) 
and we include it for generality. We perform the calculations for n = 2 k particles, where 
k E {1,..., 9}, and compute the supremum in time of the Wasserstein distance between 
subsequent solutions; cf. (3.19). We compute the Wasserstein distance by solving a linear 


programming problem based on a formulation in terms of optimal transportation. Due to the 
high computational cost (for large k), we use the following approximation 


sup W(nt 
te[o,T] 


2 fc+i 

<!h ) 


max W (/j t 
re/ 


2 fe + 1 \ II T 

>/V )='-Wk,k+ 1, 


to reduce the number of evaluations of W. Here, 


I := {jT/(N r - 1) : j = 0,..., (N r - 1)} 


(3.61) 


(3.62) 


and we take N r = 10. It should be noted, however that for the vast majority of the compu¬ 
tations, the maximum distance is observed at the final time step. 


The convergence rate for d = 1 is approximated by 

Wk+i,k- 1-2 


C Li ■= log 2 


Wk,k +1 


(3.63) 


and based on the theoretical prediction that the convergence rate is 0{n 3 ) if d = 1 , we 
expect that C\ +1 tends to the value —1. 


In Figures 3.1 3.3, results are shown for the three different values 7 E {1,2,7} respectively. 
Red graphs correspond to the scheme for 0 = 0 and blue graphs to the scheme 9 = 1. In 
these figures, the upper plots refer to computations using h = 1 for all resolutions and the 
lower plots depict computations with h varying with the number of particles used. The left 
plots show the result for density p at T = 1, as obtained with the highest resolution n = 2 9 . 
Additionally, the convergence of CW , is plotted in the right plots. 

There are several points to be mentioned about the plots. First, note that in all figures 
solutions, for h = 1 fixed and h varying, do converge to a solution by increasing the number 
of particles. The convergence is evident by the rate C^ +1 approaching its theoretical value 
— 1. Second, although convergent, solutions for h = 1 fixed and h varying are not the same 


for the same value of 7 . Third, in Figure 3.2 where 7 = 2 , the solutions obtained with the 
two schemes coincide. This effect is expected since for this value of 7 the two schemes are 
identical. On the other hand, this is not true for the the cases 7=1 and 7 = 7 . Fourth, 
interestingly enough, even though the proof only covers cases for 7 > 1 , the case 7 = 1 
converges. In the same case, it is unclear why a spike is present in the convergence graph for 
9 = 0 and varying h. Fifth, for fixed value of h all cases converge from below sharply towards 
the theoretical value Cj, +1 ~ — 1 , while for h varying with the resolution they converge from 
above. Finally, for fixed h = 1, this large value does not permit local effects to appear on 
the free boundaries of the domain. These effects are exhibited in the cases of varying h as 
discontinuities of the density profile and therefore seem to be related to problems of applying 
regularization over small h -sized regions in bounded domains. 


33 















Figure 3.3: For 7 = 7 , density p at final time T = 1, with n = 2 9 particles and convergence 
C 1 . Red and blue plots refer to the schemes for 6 = 0 and 0 = 1, respectively. Upper plots 
present results for h = 1 fixed, while lower plots to variable h = 1.5 Vq. 
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Table 1: C | +1 for the two-dimensional hydrodynamic computations 



k 

2 

3 

4 

5 

7 = 2 

9 = 

= 0 

-0.51 

-0.50 

-0.50 

-0.50 

h fixed 

9 = 

= 1 

-0.51 

-0.50 

-0.50 

-0.50 

7 = 2 

9 = 

= 0 

-0.44 

-0.47 

-0.49 

-0.44 

h varying 

9 = 

= 1 

-0.44 

-0.47 

-0.49 

-0.44 

7 = 7 

9 = 

= 0 

-0.51 

-0.50 

-0.50 

-0.50 

h fixed 

9 = 

= 1 

-0.50 

-0.50 

-0.50 

-0.50 

7 = 7 

9 = 

= 0 

-0.37 

-0.45 

-0.48 

-0.48 

h varying 

9 = 

= 1 

-0.41 

-0.43 

-0.52 

-0.51 


In two spatial dimensions, the hydrodynamic problem examined is the expansion of an initially 
square gas cloud, until time T = 1. In order to show that the results also hold for non-static 
initial conditions, a rotation described by the initial velocity field (vo, x ,vo tV ) = (—y,x) is 
applied. The same equation of state as in the one-dimensional computation is used, with 
7 £ {2,7}. The Wendland function and a constant time step of At = 10 -3 are employed. 
Note that we omit the case 7 = 1, hence do not need to ‘mimic’ Assumption |3.7| and do allow 
for bounded support in Wh- 


For d = 2, we approximate 


n 


'k +1 


:= 2 1082 


Wk+l,k+2 


Wk,k +1 


(3.64) 


Note that this definition is different from , 1 , since in d = 1 we took n of the form 2 fc , 
while in d = 2 we have n = (2 fc ) 2 , for k 6 {1, 2, 3,4, 5,6}. Here, the definition of Wk,k+i is 
modified accordingly to approximations by 4 fc and 4 fc+1 particles, respectively. The computa¬ 
tional effort for the calculation of the Wasserstein distance makes the investigation of higher n 
extremely lengthy. In the case d = 2, theory predicts that the convergence rate is 0(n -1 / 2 ), 
whence we expect that C| +1 tends to the value —1/2. In Table [lj the convergence rates 
of the two-dimensional hydrodynamic problems are shown. The theoretical value is indeed 
approached, but strong oscillations around this value appear in the case 7 = 2 with varying 
h. In Figures 3.4||3.5 particle configurations at final time T = 1 are presented for the cases 
7 = {2,7} respectively. The upper plots refer to fixed h = 1 independent of the resolution 
n (a choice in agreement with the convergence proof), while lower plots are obtained by h 
varying with the number of particles as h = 1.5i/Vo- For the plots on the left-hand side the 
scheme with 9 = 0 is used, while for the plots on the right-hand side 9 = 1 is employed. 
Similarly to the one-dimensional results, the corresponding solutions for 7 = 2 are identical 
for the schemes employing 9 = 0 or 9 = 1. On the contrary, they differ for 7 = 7 . Finally, 
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Figure 3.4: For the case 7 = 2 , particle configurations at final time T = 1 for the two- 
dimensional hydrodynamic experiment of a rotating square; on the left-hand side with 6 = 0 
and on the right-hand side with 6 = 1. The upper plots present results for h = 1 and n = 512 
particles, while the lower row refers to variable h = 1.5-v/Vo and n = 512 particles. In this 
case (7 = 2 ), the schemes for 6 = 0 and 6 = 1 are the same. 
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it should be mentioned that the instabilities of the density profile on the boundaries of the 
domain, which were observed in the one-dimensional computations, have now translated into 
the nonhomogeneous distribution of particles. 


The second numerical example considers the nonlocal force and the drag term, for which 
the numerical scheme corresponding to (3.18) does not depend on 9. Moreover, (3.18) does 
not depend on pt, hence h is only relevant if we wish to plot pt, and not for the computations 
themselves. 

The Wendland function and a constant time step of At = 10 _1 are used. For the interactions, 
we take K such that it is the gradient of the Morse potential, see e.g. [9], with parameters 
C a = 2.0, C r = 1.5, £ a = 1-0, i r = 2.0. In fact, we included a short-range regularization 
around the origin to the potential to enforce the required C*/-regularity of K. A side-effect 
is that automatically self-interactions are cancelled. Two cases for the drag coefficient are 
examined: r) = 10 and rj = 0.1, for final time T = 100. In both these cases, an equilibrium 
has been reached. Similarly to the hydrodynamical problem, (3.61) is used with n = 2 2k 
particles, where k 6 {1,2, 3, 4, 5}. Particle configurations and convergence rates are plotted 
in Figure 3.6 with the upper plots referring to 77 = 0.1 and the lower plots to rj = 10. The 
value of the convergence in this case rapidly tends to the theoretically predicted value. 


4 Concluding remarks and future directions 


Apart from the remarks already made, there are two issues that are important to point out. 
One could call them shortcomings of our approach, in the sense that these are cases to which 
our proof of convergence does not apply. The result of Theorem 3.10| does not state: 


• whether the approximations corresponding to 9 = 0 and 9 = 1, respectively, actually 
converge to the same limit solution. Our computations show that this is certainly not 
the case for h = 1 and (although the difference is smaller) neither for varying h, except 
for the trivial case 7 = 2 in which the schemes coincide. 

• whether the limit n —> 00 in any of the two cases 9 = 0 or 9 = 1 is actually like the ‘real 
physics’. To investigate this, in principle one would need to consider the limit h —> 0. 
As said before, this is beyond the scope of the current paper. 


The latter point refers to a situation in which first the limit n —* 00 is taken and afterwards 
the limit h —> 0. A more favourable approach (also from a numerical point of view) would 
be to have h depend on n in such a way that h = 0 ( 1 / yfn) as n —> 00, and hence n —> 00 
and h 


0 simultaneously. In Section |3.6 
typically done in the literature of SPH 


we anticipated this —following what is already 


and the numerical results there support the hope 
that solutions converge in the case of h varying with the number of particles. 


Nevertheless, our combined theoretical-computational results establish the convergence of 
the classical and most-used SPH scheme and also show that the corresponding equation of 
motion is a true discretized version of the equation of motion of a regularized continuous 
medium. 
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Figure 3.5: For the case 7 = 7, particle configurations at final time T = 1 for the two- 
dimensional hydrodynamic experiment of a rotating square; on the left-hand side with 6 = 0 
and on the right-hand side with 0 = 1. The upper plots present results for h = 1 and n = 512 
particles, while the lower row refers to variable h = 1.5\/Vo and n = 512 particles. 
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Figure 3.6: For the problems involving a nonlocal force term, the particle configurations at 
final time T = 100 (left plots) and the convergence rates (right plots). Drag coefficients 
77 = 0.1 (upper plot) and 77 = 10 (lower plot) are used. 
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A Completeness 


The arguments in this appendix lead to the statement of Theorem A.4 This theorem implies 
that the space C 3 defined in (3.42) is a complete metric space for every j € { 1,... ,N}. This 
result is needed to be able to apply Banach’s Fixed Point Theorem in the proof of Part [l] of 
Theorem 13.101 


Lemma A.l. Fix R > 0. Then the space V(B(R)) of probability measures on B{R ) := {.t 6 
: |x| ^ R}, endowed with the metric W := W\ , is a complete metric space. 

Proof. Since R(R) is complete, it follows from [T] , Proposition 7.1.5, that (fPi(B(R)), W\) is 
complete. Here, V\ (R(R)) is the space of probability measures with bounded first moment. 
The statement of the lemma follows from the observation that 

V 1 (B(R)) = V(B(R)). (A.l) 

The inclusion V± (R(R)) C V(B(R)) is trivial. The other inclusion follows from the fact that 
the first moment of each p 6 V(B(R)) is bounded by R. □ 

Lemma A. 2. For each T, R > 0, the space 

C([0,T];R(R(R))), (A.2) 

endowed with the metric 

sup W{m (r),^ 2 (r)), (A.3) 

re[0,T] 

is complete. 

Proof. The proof mainly follows the lines of the proof of Theorem 1.5-5 in [20j (which treats 
real-valued continuous functions). 

Let (^tn)neN denote a Cauchy sequence in C([0,T]\V(B(R))). Fix e > 0. There is a K such 
that for all m,n ^ K 

sup W(// m (r),/i n (r)) < e. (A.4) 

re[o,T] 

For any fixed t € [0,T], 

W(tl m (t),Hn{t)) C sup W{Hm{r),Hn{T)) < £ (A.5) 

re[0,T] 
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holds, so (nn(t))ne n is a Cauchy sequence in V(B(R)). It follows from Lemma A.l that 
V(B(R)) is complete and thus (li n (t)) n converges to some fit E V(B(R)). This pointwise 
limit exists for every t E [0,T], and we construct a mapping n from [0, T] to V(B(R)) by 
defining 

M*) : = h (A. 6 ) 

for all t E [0, T\. 

There is an N such that 

sup W(/r m (r),// n (r)) < e/2 (A.7) 

re[o,T] 

for all m,n ^ N (with the same e as before!). In particular, for fixed t E [0,T], 


W(limit), Unit)) < e/2 

holds for all rn. n ^ N. Thus, for each fixed t E [0, T] and for each m ^ N, 
W (ii m (t), li(t)) ^ W(ii m (t),ii n (t)) + W(ii n (t),li(t)) < e, 


(A. 8 ) 

(A.9) 


<e/2 


<e/2 


for sufficiently large n. Here we use (A. 8 ) to estimate the first term on the right-hand side. 
Due to the fact that p,(t) is defined as the pointwise limit of fi n (t), the second term can be 
made arbitrarily small by increasing n. We conclude from (A.9) that W(// m (i), fi(t)) < e for 
all m ^ N. Due to (A.7), this estimate holds with the same e and N for all t E [0, T], whence 

sup W (ii m (t), ii(t)) ^ e (A.10) 

te[o,T] 

for all m ^ N, which proves the convergence of (ii m ) to ii. 

The limit fi : [0,T] —>• V(B(R)) is continuous since it is the uniform limit of continuous 
mappings (cf. ([19] Thm. 8.3.1 for a proof for real-valued functions that can be extended 
trivially to our situation), hence ( Hn ) converges in C([0,T]; V(B(R))). □ 


Lemma A.3 (cf. [20J Theorem 1.4-7). If Y is a closed subset of a complete metric space 
(X, D), thenY is complete. 

Proof. Let (£ n ) C Y be a Cauchy sequence. Since Y C X and X is complete, there is a £ E X 
such that 

lirri d(f n . f) = 0. (A.ll) 

n—>oo 

Because (£ n ) is a sequence in Y and Y is closed, £ must be an element of Y. Thus, Y is 
complete. □ 


Theorem A.4. Define for each R > 0 

V R (W d ) := {li E V(R d ) : su PP/ u C B(R)}. (A.12) 

Fix vq E VR(R d )) and T > 0, and define 

C:={ue C([ 0 , T\;Vfi(M. d )) : i/(0) = ^o}. (A.13) 

Then the following holds: endowed with the metric 

sup W(/ii(r),^ 2 (r)), (A.14) 

re[0,T] 

the space C is a complete metric space. 
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Proof. Note that there is a one-to-one correspondence between elements of Vr{^') and ele¬ 
ments of V(B(R)). Since Lemma A.2 states that C([0,T\;V(B(R))) is complete, the same 
must hold for C([0, T]\ P^(M d )), because convergence in one of these spaces implies conver¬ 
gence in the other. We omit further details. 

Clearly, C C C([0, T\; "P^(M d )). We now show that C is closed. Let C C be a sequence 
that converges to // G C([0, T]; VR^ d ))' 


lim sup W(n(t),= 0. 
n ^°° te[o,T] 


(A.15) 


We note that 


W(/x(0), Vo) = W(/x(0), /x n (0)) ^ sup W(n(t),fj, n (t)). 

te [ o , T ] 


(A.16) 

Since the left-hand side is independent of n, while the right-hand side tends to 0 as n —» oo, 

W(/i(0), t'o) = 0 (A.17) 

must hold, so /r(0) = vq. We conclude that p € C and thus C is closed. It follows from Lemma 


A.3|that C is complete. 


a 
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